function diff = moments_SA(theta,param,aux, ind)
  
    param=param_trans_SA( theta, param);
    
	param=model_2( param,aux);
  
    if  param.K>0
        T_y = KFE_firms(param,aux);               
        ln_m=linspace(log(param.m_l),log(param.m_u),aux.n_m)';
        T_yy=T_y.*repmat((exp(ln_m)>param.m_h).*(exp(ln_m)<param.m_e),1,aux.n_m);
        T_yy=T_yy.*repmat((exp(ln_m')>param.m_h).*(exp(ln_m')<param.m_e),aux.n_m,1);

        q = q_fun(exp(ln_m),param);             
        q_p=q(2:end)-q(1:end-1);
        q_p=[q_p;0];
        h1=q_p./sum(q_p);
        Tr=expm(T_yy);
        Tr=Tr.*repmat((exp(ln_m)>param.m_h).*(exp(ln_m)<param.m_e),1,aux.n_m);
        inac=sum(Tr^12*h1);
        
    end
    
    %% Calculate moments
    u           = u_fun(param); 
    m           = linspace(param.m_l,param.m_u,aux.n)'; 
    job_to_job  = jtj_num(delta_fun(m,param),q_fun(m,param));
  
    vec = [10*(job_to_job -param.job_to_job_mom),10*(u-param.u_mom)];

    if param.K>0
       vec = [vec, inac - param.inac_mom]; 
    end     
    
    if ind == 1
        diff=vec*vec';
    elseif ind == 0
        m   = linspace( param.m_l, param.m_h, 10^3)';
        q_p = q_p_fun( m, param);
        w = q_p./length(q_p).*(param.m_h-param.m_l);
        J = J_SA( m, param);
        % Calculate implied vacancy cost. Normalize vacancy filling rate to
        % one in initial steady state (implied in chi fun)
        param.c_v = q_fun( param.m_h, param)./q_fun( param.m_u, param).*chi_fun( param).*(param.c - w'*J./q_fun( param.m_h, param));
        diff = param;
    end    
end